Article 


Water in the terrestrial planet-forming zone 


of the PDS 70 disk 


https://doi.org/10.1038/s41586-023-06317-9 


Received: 6 March 2023 


Accepted: 13 June 2023 
Published online: 24 July 2023 


Open access 


® Check for updates 


G. Perotti'™, V. Christiaens”, Th. Henning’, B. Tabone’, L. B. F. M. Waters*®, I. Kamp®, 

G. Olofsson’, S. L. Grant®, D. Gasman’, J. Bouwman’, M. Samland’, R. Franceschi’, 

E. F. van Dishoeck?”, K. Schwarz’, M. Güdel"""”?, P.-O. Lagage", T. P. Ray”, 

B. Vandenbussche”, A. Abergel*, O. Absil?, A. M. Arabhavi®, I. Argyriou®, D. Barrado”, 

A. Boccaletti"®, A. Caratti o Garatti”, V. Geers'?, A. M. Glauser’, K. Justannont”, F. Lahuis?°, 
M. Mueller®, C. Nehmé”, E. Pantin”, S. Scheithauer', C. Waelkens®, R. Guadarrama”, H. Jang’, 
J. Kanwar®7"22, M. Morales-Calderon", N. Pawellek", D. Rodgers-Lee”, J. Schreiber’, 

L. Colina”, T. R. Greve”, G. Östlin” & G. Wright® 


Terrestrial and sub-Neptune planets are expected to form in the inner (less than 10 Au) 
regions of protoplanetary disks!. Water plays a key role in their formation? *, although 
itis yet unclear whether water molecules are formed in situ or transported from the 
outer disk*ć. So far Spitzer Space Telescope observations have only provided water 
luminosity upper limits for dust-depleted inner disks’, similar to PDS 70, the first 
system with direct confirmation of protoplanet presence®’. Here we report JWST 
observations of PDS 70, a benchmark target to search for water ina disk hosting a 


large (approximately 54 AU) planet-carved gap separating an inner and outer dis 


Ke” 


Our findings show water in the inner disk of PDS 70. This implies that potential 
terrestrial planets forming therein have access to a water reservoir. The column 
densities of water vapour suggest in-situ formation via a reaction sequence involving 
O, H, and/or OH, and survival through water self-shielding>. This is also supported 

by the presence of CO, emission, another molecule sensitive to ultraviolet 
photodissociation. Dust shielding, and replenishment of both gas and small dust 
from the outer disk, may also play a role in sustaining the water reservoir”. Our 
observations also reveal a strong variability of the mid-infrared spectral energy 
distribution, pointing to a change of inner disk geometry. 


Observations of PDS 70 were taken with the JWST Mid-InfraRed Instru- 
ment (MIRI)? Medium Resolution Spectrometer” (MRS; spectral 
resolving power R = 1,600-3,400) as part of the guaranteed time MIRI 
mid-INfrared Disk Survey (MINDS; see Methods and Extended Data 
Fig. 1). The complete spectrum of PDS 70 shows several distinct traits 
(Fig. 1), which stands out with respect to other T Tauri disks”. 

Asignificant flux offset—up to a factor of 1.5 at wavelengths beyond 
18 pm—is found between the MIRI and the archival Spitzer Infra- 

Red Spectrograph (IRS) low-resolution (R = 60-100) spectra recorded 
with 15 years and one day time difference. This discrepancy is too large 
to be explained by calibration uncertainties; the absolute uncertainty 
for both IRS and MIRI is approximately 5% for the 4.9-22.5 um range”. 
Similarly, the difference in aperture size of the two spectrographs 


cannot account for such an offset. Hence, with the current MIRI 
data reduction, time variability is the most likely explanation for the 
observed flux differences. 

Variability in the mid-infrared observed with Spitzer-IRS has 
been mainly attributed to short-wavelength stellar irradiation or to 
dynamical changes in the inner disk geometry due to the presence 
of planets’®”°. In the case of PDS 70, stellar irradiation is excluded 
as it would cause an overall increase or decrease in flux, contrary to 
what is observed with the Wide-field Infrared Survey Explorer (WISE) 
time-series observations (Extended Data Fig. 6). PDS 70 is known to 
be ina late stage of accretion—with an estimated waning accretion 
rate™ of approximately 10™ M, yr‘ making it unlikely to explain the 
significant flux difference. Changes in the scale height of the inner 
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Fig. 1|JWST-MIRI MRS spectrum of PDS 70. The spectrum is a composite of 
the colour-coded short, medium and long sub-bands of the four MIRIMRS 
Integral Field Units (IFUs)”. The Spitzer-IRS spectrum of PDS 70 is also shown 
in grey. The major dust features are labelled. The spectrum is dominated by 
exceptionally prominent silicate emission at 10 and 18 um and it clearly shows 
anumber of crystalline dust features. The much higher sensitivity and spectral 


disk wall emitting at shorter wavelengths (approximately 2-8 um) 
can be responsible for shadowing the disk material located further 
out, resulting in less emission at wavelengths beyond 18-20 um (also 
referred to as ‘seesaw-like’ variability”). However, for PDS 70 a complete 
seesaw-like profile is not observed as there is no corresponding increase 
in flux at the shorter wavelengths. Time variability is also supported 
by the above-mentioned WISE observations, which indicate that such 
variability occurs on short timescales (of at most 1 yr) and that it may 
indeed be attributed to occulting material located close to the star 
(approximately 1 AU). 

The MIRI spectrum of PDS 70 clearly shows the presence of silicate 
dust grains that have undergone significant thermal processing (Fig. 1). 
We attribute the crystalline dust features to enstatite at 9.40 um and 
forsterite at 11.30 and 16.40 um. The observed dust continuum is well 
reproduced with a three-component disk model, with a 400-600 K 
surface layer accounting for the bulk of the observed emission (Fig. 2). 

The MIRI spectrum also reveals a wealth of water lines, particularly 
in the 7 pm spectral window (Fig. 3). This indicates the presence of a 
water reservoir in the terrestrial region of a disk already hosting two or 
more protoplanets. As such, it also provides important clues to theories 
on the origin of water during terrestrial planet formation in the solar 
system?**, We focus on the ro-vibrational transitions of the bending 
mode of para- and ortho-water in the 7 um region where the brightest 
lines are observed and contamination by the stellar atmosphere is 
negligible (Extended Data Fig. 3). This includes strong water blends 
dominated by lines with upper energy level E, ~ 2,400-3,200 K. Weaker 
lines are also detected at the 1 mJy level, some of them corresponding 
to more excited levels up to E, ~ 4,300 K. 

Further insight into the origin of water emission is obtained from 
zero-dimensional slab modelling, which has also been used to inter- 
pret Spitzer spectroscopic data”. The synthetic spectrum of water is 
calculated froma plane-parallel slab model, where the level populations 
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resolution of MIRI MRS compared with Spitzer-IRS allows us to detect an inner 
disk gas reservoir by showing weak emission of water vapour and carbon dioxide 
as well as two molecular hydrogen lines. The insets show the ro-vibrational and 
rotational transitions of ortho- and para-H,O, the molecular hydrogen H,S(1) 
and S(5) rotational lines, and the v; bending mode of CO). 


are in local thermodynamic equilibrium (LTE) at a single excitation 
temperature T. The other fitting parameters are the line-of-sight 
column density Nwithin an effective emitting area nR’ given by its radius 
R, and the intrinsic line broadening assumed to be g = 2 km s” (ref. 25). 
Note that R does not need to correspond toa disk radius, but could also 
represent an annulus with the same area or an emission spot breaking 
the axisymmetry. The best-fit model is then obtained by minimizing 
the reduced x” between measured and model line fluxes over the indi- 
vidual spectral window around each H,O line (Extended Data Fig. 4). 

The observed H,O spectrum in the 6.78-7.36 um spectral region is 
best fitted witha slab of gas at T= 600 K, with an emitting area of radius 
R=0.047 au and a column density of N=1.4 x 10 cm”. The tempera- 
ture is mostly determined by the ratio between the lines of different 
E,, for example, the series of lines in the 7.3 um region. The column 
density is set by the ratio between the weaker lines and indicates that 
the brightest lines are optically thick. The emitting area is constrained 
by matching the fluxes of the optically thick lines, and points towarda 
compact emission region. This is further supported by the fact that the 
detected lines are broad (AA = 0.01-0.05 um); if the line broadening is 
caused by the gas kinematics, the full-width half-maximum (FWHM) 
of the line would be about 100 kms“, corresponding to a Keplerian 
radius of 0.1 AU, consistent with the emitting area deduced from our 
fit after correction for disk inclination i= 51.7 + 0.1 (ref. 11). Interest- 
ingly, we find that our best-fit LTE model of the water emission in the 
7 um region reproduces water rotational lines at 15 um reasonably well, 
suggesting that all water emission inthe MIRI spectral range originates 
from inside approximately 0.05 AU under LTE conditions (Extended 
Data Fig. 5). 

Besides water in the 7, 15 and 17 pm regions, other species have been 
identified but the analysis is postponed to a future study. The fun- 
damental Q-branch of CO, corresponding to the v; bending mode is 
detected at 14.96 um (Extended Data Fig. 5). Interestingly, the width 
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Fig. 2 | Dust continuum fit to the MIRI spectrum of PDS 70. a, The disk model 
has three spectral components: an inner rim, an optically thick midplane disk 
layer and an optically thin warm disk surface layer. The stellar photospheric 
emission is represented by a stellar atmosphere model (see Methods for 
further details). The surface layer component dominates the MIRI spectrumin 
the 8-22.5 um wavelength range. Its temperature is constrained to be between 


Continuum-subtracted flux (mJy) 


400 and 600K. The silicate emission at 8-12 um is consistent with a population 
ofoptically thin dust grains with typical sizes of 0.1-2 um. A significant 
contribution from an optically thick dust component is excluded because of 
the high silicate peak/continuumratio of approximately 4 (ref. 31). b, The 
residuals on the dust continuum fit. 
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Fig. 3 | Continuum-subtracted spectrum showing H,O emission in the 7 pm region and the best-fit LTE slab model. The best-fit model (blue) has T= 600 K, 
N(H,O)= 1.4 x 108cm? and R = 0.047 au. The molecular hydrogen H, S(5) line is labelled on top of the spectrum. 
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Fig. 4| Comparison between water luminosity and mid-infrared spectral 
index (3-39) for a sample’ of protoplanetary disks. n,;_;, isa diagnostic of 
the presence and size of inner disk dust cavities: 0.9 < n,3_39 < 2.2 corresponds 
to disks with large gaps and/or cavities**”’. Black dots represent disks with 
mid-infrared water detections. Disks for which only upper limits were obtained 
are shown as gold arrows. The grey shaded area highlights the location of 
dust-depleted inner disks with PDS 70 shown as a blue dot. The water line flux 
used to compute the water luminosity of PDS 70 is calculated as describedina 
previous work’, The Spitzer spectrum is used to estimate n,;_4, for PDS 70 to be 
consistent with the other targets. Spitzer-IRS obtained only water luminosity 
upper limits for disks characterized by n,3-3) greater than 0.9. Below 10 um, IRS 


of this feature is sensitive to the temperature and indicates cooler gas 
at T~ 200 K inthe optically thin regime. The pure rotational molecular 
hydrogen H,S(5) and H,S(1) lines are detected at 6.91 and 17 um (Fig. 1). 
We note that the H, S(2), S(3) and S(4) lines coincide with the broad 
silicate emission feature and thus establishing their presence needs 
to await an in-depth analysis of this dust feature. 

Spitzer-IRS observations detected water in approximately 50% of 
dust-rich inner disks around T Tauri stars”°, but obtained only upper 
limits for disks with large inner dust gaps or cavities defined by a 
mid-infrared spectral index 13-3, > 0.9, where nis the slope of the spec- 
trum between 13 and 30 um (Fig. 4)’. The detection of water vapour in 
the PDS 70 MIRI spectrum demonstrates that PDS 70 has maintained 
tosome degree the physical and chemical conditions of dust-rich inner 
disks in its terrestrial planet-forming zone despite the presence of a 
notably large gap (see ‘Origin of water in PDS 70’ in Methods). Our LTE 
slab model only provides a first quantitative analysis of the H,O emission. 
Non-LTE effects could lead to subthermal line emission, which would 
make our estimated emitting area a lower limit. In T Tauri disks with 
strong radial temperature gradients, the water lines are expected to orig- 
inate from different regions of the disk depending on their upper energy 
level and Einstein-A coefficients”. Detailed modelling using a realistic 
disk structure and including non-LTE effects such as infrared radiative 
pumping is needed in the future to further constrain the distribution 
of H,O across the inner disk. However, this first analysis already proofs 
that the inner disk of PDS 70 is rich in water and the inferred slab model 
parameters are roughly consistent with a detailed thermo-chemical 
model (ref. 28 and B. Portilla-Revelo, personal communication). 

The luminosity of the 17 um water lines is two orders of magnitude 
weaker for PDS 70 than DoAr 44 (refs. 7,29). DoAr 44 is a system with 
similar properties to PDS 70, but characterized by 1,3_3. < 0.9 (Fig. 4). 
This result points to a colder water reservoir in PDS 70, and is con- 
sistent with the lower luminosity and lower accretion rate of PDS 70 
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provided only a spectral resolution of R ~100, preventing a comprehensive 
view of water in the innermost regions**. DoAr 44 is a system schematically 
similar to PDS 70 (ref. 29). The two stars have comparable age and spectral 
types K3 and K7, respectively; DoAr 44 has a higher mass accretion rate of 

Maco *™ 10° Mo yr! (ref. 30). The cavity size of DoAr 44 is 34 Au (refs. 7,29), which 
is smaller than that of PDS 70 (approximately 54 Au (ref. 11)). Both systems have 
small inner disks based on VLTI-GRAVITY*, VLT-SPHERE and ALMA data®", and 
in both systems the water emission is contained to within 1 Au (ref. 29). The 
water luminosity of PDS 70 is two orders of magnitude weaker than that of 
DoAr 44, pointing toa colder water reservoir in PDS 70. This is consistent with 
the lower luminosity and lower accretion rate of PDS 70. 


(refs. 21,30). This work opens anew window onthe origin of water in pro- 
toplanetary disks by showing that MIRI MRS can now detect very weak 
(s5 mJy) water lines in the innermost regions of disks with large gaps, 
and hence that the presence of water in the terrestrial planet-forming 
zone of dust-depleted inner disks is not as rare as previously thought. 


Online content 


Any methods, additional references, Nature Portfolio reporting summa- 
ries, source data, extended data, supplementary information, acknowl- 
edgements, peer review information; details of author contributions 
and competing interests; and statements of data and code availability 
are available at https://doi.org/10.1038/s41586-023-06317-9. 


Mulders, G. D., Pascucci, l., Apai, D. & Ciesla, F. J. The Exoplanet Population Observation 
Simulator. |. The inner edges of planetary systems. Astron. J. 156, 24 (2018). 

2. Ciesla, F. J. & Cuzzi, J. N. The evolution of the water distribution in a viscous 
protoplanetary disk. Icarus 181, 178-204 (2006). 
Eistrup, C. & Henning, T. Chemical evolution in ices on drifting, planet-forming pebbles. 
Astron. Astrophys. 667, A60 (2022). 

4. Krijt, S. et al. Chemical habitability: supply and retention of life's essential elements 
during planet formation. Preprint https://arxiv.org/abs/2203.10056 (2022). 
Bethell, T. & Bergin, E. Formation and survival of water vapor in the terrestrial 
planet-forming region. Science 326, 142-153 (2009). 

6. Glassgold, A. E. et al. Formation of water in the warm atmospheres of protoplanetary 
disks. Astrophys. J. 701, 1675-1677 (2009). 
Banzatti, A. et al. Hints for icy pebble migration feeding an oxygen-rich chemistry in the 
inner planet-forming region of disks. Astrophys. J. 903, 124 (2020). 

8. Keppler, M. et al. Discovery of a planetary-mass companion within the gap of the 
transition disk around PDS 70. Astron. Astrophys. 617, A44 (2018). 
Haffert, S. Y. et al. Two accreting protoplanets around the young star PDS 70. Nat. Astron. 
3, 749-754 (2019). 

10. Long, Z. C. et al. Differences in the gas and dust distribution in the transitional disk of a 
Sun-like young star, PDS 70. Astron. Astrophys. 858, 112 (2018). 
Keppler, M. et al. Highly structured disk around the planet host PDS 70 revealed by 
high-angular resolution observations with ALMA. Astron. Astrophys. 625, A118 (2019). 

12. Benisty, M. et al. A circumplanetary disk around PDS70c. Astrophys. J. Lett. 916, L2 (2021). 


13. 


14. 


15. 


16. 


17. 


18. 


19. 


20. 


21. 


22. 


23. 


24. 


25. 


26. 


27. 


28. 


Rieke, G. H. et al. The Mid-Infrared Instrument for the James Webb Space Telescope, I: 
introduction. Publ. Astron. Soc. Pac. 127, 584 (2015). 

Wright, G. S. et al. The Mid-Infrared Instrument for the James Webb Space Telescope, Il: 
design and build. Publ. Astron. Soc. Pac. 127, 595 (2015). 

Wells, M. et al. The Mid-Infrared Instrument for the James Webb Space Telescope, VI: the 
medium resolution spectrometer. Publ. Astron. Soc. Pac. 127, 646 (2015). 
Kessler-Silacci, J. et al. c2d Spitzer IRS spectra of disks around T Tauri stars. |. Silicate 
emission and grain growth. Astrophys. J. 639, 275-291 (2006). 

Furlan, E. et al. A survey and analysis of Spitzer infrared spectrograph spectra of T Tauri 
Stars in Taurus. Astrophys. J. 165, 568-605 (2006). 

Argyriou, l. et al. JWST MIRI flight performance: the medium-resolution spectrometer. 
Preprint at https://arxiv.org/abs/2303.13469 (2023). 

Muzerolle, J. et al. Evidence for dynamical changes in a transitional protoplanetary disk 
with mid-infrared variability. Astrophys. J. Lett. 704, L15-L19 (2009). 

Espaillat, C. et al. A Spitzer IRS study of infrared variability in transitional and 
pre-transitional disks around T Tauri stars. Astrophys. J. 728, 49 (2011). 

Manara, C. F. et al. Constraining disk evolution prescriptions of planet population 
synthesis models with observed disk masses and accretion rates. Astron. Astrophys. 
631, L2 (2019). 


Skinner, S. L. & Audard, M. HST UV spectroscopy of the planet-hosting T Tauri star PDS 70. 


Astrophys. J. 938, 134 (2022). 
Stimpfl, H. An angstrém-sized window on the origin of water in the inner Solar System: 


Atomistic simulation of adsorption of water on olivine. J. Cryst. Growth 294, 83-95 (2006). 


Genda, H. & Ikoma, M. Origin of the ocean on the Earth: early evolution of water D/H ina 
hydrogen-rich atmosphere. Icarus 194, 42-52 (2008). 
Salyk, C. et al. A Spitzer survey of mid-infrared molecular emission from protoplanetary 


disks. II. Correlations and local thermal equilibrium models. Astrophys. J. 731, 130 (2011). 


Pontoppidan, K. M. et al. A Spitzer survey of mid-infrared molecular emission from 
protoplanetary disks. |. Detection rates. Astrophys. J. 720, 887-903 (2010). 

Blevins, S. M. et al. Measurements of water surface snow lines in classical protoplanetary 
disks. Astrophys. J. 818, 22 (2016). 

Portilla-Revelo B. et al. Constraining the gas distribution in the PDS 70 disk as a method to 
assess the effect of planet-disk interactions. Preprint at https://arxiv.org/abs/2306.16850 
(2023). 


29. Salyk, C. et al. Detection of water vapor in the terrestrial planet forming region of a 
transition disk. Astrophys. J. Lett. 810, L24 (2015). 

30. Manara, C. F. et al. Gas content of transitional disks: a VLT/X-Shooter study of accretion 
and winds. Astron. Astrophys. 568, A18 (2014). 

31. Oliveira, I. et al. A Spitzer Survey of protoplanetary disk dust in the Young Serpens Cloud: 

how do dust characteristics evolve with time? Astrophys. J. 714, 778-798 (2010). 

32. Brown, J. M. et al. Cold disks: Spitzer spectroscopy of disks around young stars with large 

gaps. Astrophys. J. Lett. 664, L107-L110 (2007). 

33. Furlan, E. et al. Disk evolution in the three nearby star-forming regions of Taurus, 

Chamaeleon, and Ophiuchus. Astrophys. J. 703, 1964-1983 (2009). 

34. Banzatti, A. et al. The kinematics and excitation of infrared water vapor emission from 

planet-forming disks: results from spectrally-resolved surveys and guidelines for JWST 

spectra. Astron. J. 165, 72 (2023). 

35. Bouvier, J. et al. Investigating the magnetospheric accretion process in the young 

pre-transitional disk system DoAr 44 (V2062 Oph). A multiwavelength interferometric, 

spectropolarimetric, and photometric observing campaign. Astron. Astrophys. 643, A99 

(2020). 


Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in 
published maps and institutional affiliations. 


Open Access This article is licensed under a Creative Commons Attribution 

ny 4.0 International License, which permits use, sharing, adaptation, distribution 

and reproduction in any medium or format, as long as you give appropriate 

credit to the original author(s) and the source, provide a link to the Creative Commons licence, 
and indicate if changes were made. The images or other third party material in this article are 
included in the article's Creative Commons licence, unless indicated otherwise in a credit line 
to the material. If material is not included in the article's Creative Commons licence and your 
intended use is not permitted by statutory regulation or exceeds the permitted use, you will 
need to obtain permission directly from the copyright holder. To view a copy of this licence, 
visit http://creativecommons.org/licenses/by/4.0/. 


© The Author(s) 2023 


Nature | www.nature.com | 5 


Article 


Methods 


PDS 70 system 

PDS 70 (V1032 Cen) is a K7-type star in the Upper Centaurus-Lupus 
subgroup (d= 113.4 + 0.5 pc (ref. 36)) in a late stage of accretion” 
with an estimated age of 5.4 +1.0 Myr (ref. 37). The disk around PDS 
70 (refs. 38-40) hosts two actively accreting protoplanets: PDS 70 b 
and PDS 70 c, which reside in an approximately 54 AU annular gap 
between an inner and outer disk®’. The presence of an inner dusty 
disk in the PDS 70 system has been inferred from both near-infrared 
scattered light and ALMA images” ”. The 855 um dust continuum 
emission from the innermost disk regions is confined within the orbit 
of PDS 70 b (approximately 22 au; Extended Data Fig. 1), putting an 
upper limit to the inner disk radial extent of approximately 18 AU 
(ref. 12). A population of small dust grains may be responsible for the 
observed inner disk emission although the current low dust mass esti- 
mates could support the simultaneous presence of small and large 
dust grains“. 


Observations and data reduction 

The PDS 70 disk (CD-40-8434) was observed with MIRI?” on1 August 
2022 as part of the Guaranteed Time Observation (GTO) programme 
1282 (PI: Th. Henning) with number 66. The disk was observed in 
FASTR1 readout mode with a four-point dither pattern in the nega- 
tive direction for a total on-source exposure time of 4,132 s. The MRS” 
mode was used, which has four IFUs. Each IFU (referred to as chan- 
nel) covers a different wavelength range and splits the field of view 
into spatial slices. Calibration and processing of IFU observations 
produces three-dimensional spectral cubes. The latter are used to 
extract a final spectrum covering the MIRI 4.9-22.5 um range and isa 
composite of the four IFUs: channel 1 (4.9-7.65 um; R = 3,400), channel 
2 (7.51-11.71 um; R = 3,000), channel 3 (11.55-18.02 um; R = 2,400) and 
channel 4 (17.71-22.5 um; R = 1,600). Each channel is in turn composed 
of three sub-bands: SHORT (A), MEDIUM (B) and LONG (C) leading to 
a total of 12 wavelength bands. 

We processed the PDS 70 data using a hybrid data reduction pipe- 
line made from the combination of the JWST Science Calibration 
pipeline” (v.1.8.4) stages 1 to 3, with dedicated routines based on the 
Vortex Image Processing (VIP) package**“* for bad pixel correction, 
background subtraction and removal of spikes affecting the final 
spectrum. Specifically, data reduction proceeded as follows: (1) the 
class Detectorl of the WST pipeline was used to process uncalibrated 
raw data files using Calibration Reference Data System (CRDS) con- 
text jwst_1019.pmap and default parameters; (2) apart from pixels 
flagged inthe Data Quality (DQ) extension, we identified additional bad 
pixels with both an iterative sigma clipping algorithm and througha 
cross-shaped match filter, and corrected them using a two-dimensional 
Gaussian kernel; (3) Spec2 was then used with default parameters, but 
the background subtraction was skipped, and dedicated reference 
files” for photometric and fringe flat calibrations were adopted; (4) 
as no dedicated background observation was taken, we leveraged the 
four-point dither pattern to obtain a first guess on the background 
map, then refined it using a median filter, which both smoothed the 
background estimate and removed residual star signals from it; (5) 
Spec3 was then run with default parameters, apart from the mas- 
ter_background and outlier_detection steps that were turned off, 
in the latter case to avoid spurious spectral features resulting from 
under-sampling of the point spread function; (6) we recentred the 
spectral cubes by applying the shifts maximizing the cross-correlation 
between cube frames, and found the location of the point spread func- 
tion centroid with a two-dimensional Gaussian fit on the median image 
of each aligned cube; (7) spectra were then extracted with aperture 
photometry in 2.5-FWHM apertures centred onthe centroid location 
(with the FWHM equal to approximately 1.22A/D with the telescope 
diameter D equal to 6.5 m), corrected for both aperture size using 


correction factors’, and spikes affecting individual spaxels included 
inthe aperture; (8) spectra were finally corrected for residual fringes 
at the spectrum level and the bands were stitched together based on 
the level of the shorter wavelength bands (these rescaling factors were 
systematically within 3% of the photometric solution). Spurious data 
reduction artefacts were masked at 5.12, 5.90, 7.45 and 7.50 um. The 
uncertainty associated with each photometric measurement consid- 
ers both Poisson and background noise, combined in quadrature. The 
former is an output from theJWST pipeline, whereas for the latter we 
propagated our background estimate obtained in step (4) through 
Spec3, and considered the standard deviation of the fluxes inferred 
in independent 2.5-FWHM apertures as a proxy for the background 
noise uncertainty. The final relative uncertainties range from approxi- 
mately 0.1% to approximately 1.1% with respect to the continuum at 
the shortest and longest wavelengths considered in this work (4.9 and 
22.5 um, respectively). 


Local continuum fit 

Extended Data Fig. 2 shows the local baseline fit for the 7 um region. 
The continuum level is determined by selecting line-free regions and 
adopting a cubic spline interpolation (scipy.interpolate.interpld). 
This continuum is then subtracted from the original data to produce 
the spectrum shown in Fig. 3. 


Correction for the photospheric emission 

The observed near-infrared colour index,/—K,=1.01(2MASS), indicates 
asmall colour excess, E(J - K) = 0.16, which could be due to either inter- 
stellar extinction or a true excess in the K, band, or a combination of the 
two. By assuming that the brightness in the K, band is essentially due to 
photospheric emission we can, by using a model atmosphere provided 
by P. Hausschildt (personal communication, 2023), extrapolate the con- 
tribution into the mid-infrared spectral region. The parameters used 
for the model atmosphere of the PDS 70 K7-star are an effective tem- 
perature 7,,,= 4,000 K, surface gravity log(g) = 4.5 and solar metallicity. 
At5 pm the photospheric contribution amounts to 56 mJy, that is, 44% 
of the observed flux density. At longer wavelengths, inthe 7 um region 
in which the water emission is detected, the photosphere amounts to 
one-third of the observed flux density and the subtraction of the pho- 
tospheric contribution marginally alters the continuum-subtracted 
spectrum (Extended Data Fig. 3). 


Slab models fits 

The molecular lines are analysed using a slab approach that takes into 
account optical depth effects. The level populations are assumed to 
be in LTE and the line profile function to be Gaussian with a FWHM of 
AV=4.7 kms” (o =2 km s)”. The line emission is assumed to originate 
from aslab of gas with a temperature 7 and a line-of-sight column den- 
sity of N. Under these assumptions and neglecting mutual line opacity 
overlap, the frequency-integrated intensity of a line is computed as 
follows: 


AY 


~ 2,/In2A4 


where B, (T)is the Planck function, Ay is the rest wavelength of the line 
and Tois the optical depth at the line centre vo, with: 
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In this equation, x, and x, denote the level population of the lower 

and upper states, g; and g, their respective statistical weights and A,, 

the spontaneous downward rate of the transition. The line intensity is 


then converted into integrated flux F, assuming an effective emitting 
area of nR? anda distance to the source das: 


anf — exp(- we” fy (1) 


2 
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We note that neglecting mutual line overlap for H,O when calculating 
the line intensity is a valid approximation for N(H,O) < 10° cm“ and 
significantly reduces the computational time*ć. Finally, the spectrum 
is convolved and sampled in the same way as the observed spectrum” 
and all lines are then summed to prepare a total synthetic spectrum. 
The molecular data, that is, line positions, Einstein A coefficients and 
statistical weights stem froma previous work*’. 


Fitting procedure for H,O vapour lines 

The LTE slab model described above is then used to fit the HO lines 
in the 6.78-7.36 um region following ay’ method. First, an extended 
grid of models is computed varying the total column density from 10” 
to 10” cm” in steps of 0.17 in log,, space and the temperature from 
100-1,400 K in steps of 50 K. We further assume an ortho-to-para ratio 
of 3. For each set of free parameters (N, T, R), a synthetic spectrum is 
calculated at the spectral resolving power R = 2,000 and rebinned to 
the spectral sampling of the observed spectrum using the slabspec 
python code”. The adopted spectral resolution is lower that the nomi- 
nal MIRI MRS spectral resolution in channel 1 (ref. 15); it was selected 
to account for the observed line broadening. This spectrum is fur- 
ther used to compute the x’ value on a spectral channel basis. Spe- 
cific spectral windows are chosen to avoid contamination by other 
gas features. This includes all spectral channels falling within 0.02 um 
(1,000 km s”) of any hydrogen recombination line and a 0.01 um wide 
spectral windowat the position of the S(5) line of H, at 6.91 um. To mit- 
igate the errors induced by the continuum subtraction procedure, 
we also include only spectral elements falling within 0.004 um ofa 
water line. For each value of (N, T), the x? is then minimized by vary- 
ing the emitting size nR’. The resulting x? map is shown in Extended 
Data Fig. 4 together with the best-fit emitting radius. The confidence 
intervals are estimated following a previous work® and using a repre- 
sentative noise level of ø = 0.15 mJy. We note that, owing to the large 
number of lines in this crowded region, there is little space to deter- 
mine the noise on the continuum. Therefore, we estimate the noise 
level between 7.72 um and 7.73 um to avoid contamination by H,O and 
hydrogen lines. 


Origin of water in PDS 70 

At the typical densities of inner disk regions (n,,> 10° cm) the chem- 
istry can rapidly reach steady-state conditions and water vapour can 
form froma simple reaction sequence involving O, H, and OH. Water 
and OH absorb efficiently in the ultraviolet (that is, water and OH 
shielding), ensuring the survival of water molecules even in regions 
of reduced dust opacity’. This mechanism by itself is able to account 
for the column densities of water vapour detected in this work and is 
supported by the presence of CO, emission. Small grains in the inner 
disk provide additional ultraviolet shielding. One question that natu- 
rally arises is whether the water vapour in PDS 70 originated before the 
formation of the giant protoplanets within the gap or whether there 
is acontinuous supply of gas from the outer to the inner disk regions. 
ALMA high spatial resolution CO observations reveal the presence of 
gas inside the gap"*". Observations and models find the gap to be gas 
depleted" (two to three orders of magnitude assuming an r* surface 
density profile) and dust depleted”, but not empty. One possibility 
could be that a population of water-containing dust particlesis able to 
filter through the orbits of PDS 70 b and PDS 70c, enriching the inner 
disk reservoir”. Experimental evidence indicates that water chemically 
bound to complex silicates can be preserved to temperatures up to 
400-500 K (refs. 53,54) and thus survive in the regions probed by our 
observations inside the water snowline. We note that some degree of 
dust filtering is expected with gas replenishment, as small dust particles 


can couple tothe gas. Therefore, a replenishment of both gas and dust 
from the outer disk to sustain the water reservoir and hence the PDS 
70 accretion rate is possible. 


Fitting procedure for the dust continuum 

The 4.9-22.5 um dust continuum is analysed using a two-layer disk 
model for the dust emission”. This model was successfully applied to 
Spitzer-IRS spectra of planet-forming disks*; we follow the same mod- 
elling approach here. We rebin the spectrum by averaging 15 spectral 
points and assign errors oto the rebinned spectral points assuming a 
normal error distribution with equal weights for each individual spec- 
tral element. The stellar photospheric emission is represented by a stel- 
lar atmosphere model fitted to optical and near-infrared photometry. 
The disk model has three spectral components: (1) a hot inner disk 
Fim (2) an optically thick midplane disk layer Fn and (3) an optically 
thin warm disk surface layer F,,,,. The dust grains representing the disk 
components are assumed to have power-law temperature distributions, 
and each is characterized by a minimum anda maximum temperature 
Tam. The disk surface layer is assumed to consist of a number of dust 
species iwith different chemical compositions and witha fixed number 
of grain sizes /, all emitting at the same temperatures. The total disk 
flux can then be written as: 


Fy = Fy timt Fv,mp * Ay, sur (4) 
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and B,(7) is the Planck function, qam is the power-law exponent for the 
temperature gradient in the disk surface layer, x; jare the opacities in 
cm? g” of dust species i with grain size j, dis the distance to the star 
and D, jare normalization factors”. We use three grain compositions 
(with SiO,, SiO, and SiO, stochiometry)*” “ and both amorphous and 
crystalline lattice structures to capture the rich spectral structure evi- 
dent in the MIRI data. The choice of this set of compositions is based 
on previous analyses”, which showed that this set of materials is able 
to capture most spectral variations in planet-forming disks observed 
with Spitzer-IRS. We use either two or three grain sizes (that is, 0.1, 2 
and 5 um) for each of the dust species. In total, the model has 23 fitting 
parameters. We use the MultiNest Bayesian fitting algorithm” and the 
PyMultiNest package® to find the best-fit parameters. The resulting fit 
and the separate spectral components (star, inner rim, midplane and 
surface layer) are shown in Fig. 2. 


WISE time-series observations 

Extended Data Fig. 6 reports WISE time-series observations of PDS 70. 
Observations were executed on 2-3 February and 6 February 2010, and 
on1-2 August 2010. We note that the source is highly variable and that 
WISE 4 (W4; 25 um) is anticorrelated with WISE 1 (W1; 3.4 ym) and WISE 2 
(W2; 4.6 um). Such variability may be ‘seesaw’-like”, for which changes 
inthe scale height of the inner disk wall shadow the disk material located 
further out. However, acomplete ‘seesaw’ profile is not observed, as at 
wavelengths shorter than 8 pm the MIRI spectrum lies above the IRS 
spectrum (Fig. 1). This is not surprising as the wavelength of the ‘pivot’ 
point (that is, the wavelength at which a shift in emission is observed) 
is dependent on the location of the occulting material with respect to 
the star, the stellar luminosity and the inclination of the system, with 
highly inclined systems showing a more complete ‘seesaw’ than more 
face-on systems such as PDS 70 (i= 51.7 + 0. 1 (ref. 11)). Interestingly, 
WISE 3 (W3;12 pm) is not anticorrelated with WISE 1 and WISE 2 because 
of the dominant 10 um silicate band, which indeed shows a minor offset 
compared with the longer wavelengths. This indicates that the material 
contributing to the 10 um emission is not shadowed. This behaviour is 


Article 


seen if the emission arises from warmer dust closer to the star than the 
occulting material or further above the disk midplane. 

We also note that the difference in aperture size of Spitzer-IRS and 
MIRI MRS cannot explain the observed variability. The Spitzer-IRS 
low-resolution spectrograph has a slit width of 3.6” for wavelengths 
shorter than 14 um and 10.2” for wavelengths longer than 14 pm. 
Although the maximum aperture of Spitzer-IRS at longer wavelengths 
is larger than that of MIRI MRS, this is not the case for the shorter wave- 
lengths, for which the slit widths are similar for both observatories 
(approximately 3.6” versus 4.0”). However, a flux offset is also observed 
inthis spectral region. In addition, in the case when the long wavelength 
excess would arise from an extended component, ajump in flux level 
at 14 um—where the aperture size changes—would be present in the 
IRS data, but it is absent. 


Data availability 


The original data analysed in this work are part of the GTO programme 
1282 (PI: Th. Henning) with number 66 and will become public on 
2 August 2023 onthe MAST database (https://mast.stsci.edu). The por- 
tion of the spectrum presented in Fig. 3 is available on Zenodo at https:// 
zenodo.org/record/7991022. The spectroscopic data for water can 
be downloaded from the HITRAN database (https://hitran.org). The 
Spitzer-IRS spectrum plotted in Fig. Lis part of the Spitzer-IRS GTO pro- 
gramme 40679 (PI:G. Rieke). The spectrum was extracted and calibrated 
using private codes and is available on Zenodo at https://zenodo.org/ 
record/7991022. The optical constants of the dust species considered 
in the fitting procedure for the dust continuum can be downloaded 
from the HJPDOC database (https://www2.mpia-hd.mpg.de/HJPDOC). 


Code availability 


The slab model used in this work is a private code developed by B.T. and 
collaborators. It can be obtained from B.T. upon request. The synthetic 
spectra presented in this work can be reproduced using the slabspec 
code, which can be found at https://doi.org/10.5281/zenodo0.4037306. 
The fitting procedure for the dust continuum uses the publicly available 
MultiNest Bayesian fitting algorithm (https://github.com/Johannes 
Buchner/MultiNest) and the PyMultiNest package (https://github.com/ 
JohannesBuchner/PyMultiNest). Figures were made with Matplotlib 
v.3.5.1. under the Matplotlib license at https://matplotlib.org/. 
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Extended Data Fig. 2 | Local continuum fit used in the spectrum presented in Fig. 3. a, The selected continuum points are displayed as red dots and the 
interpolated continuum is shown as a red line. b, The continuum-subtracted spectrum. 
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Extended Data Fig. 3| Correction for the photospheric emission. 
a, Comparison between the MIRI-MRS spectrum (black) and the spectrum 
corrected for the stellar photosphere (orange). Both spectra are continuum 


subtracted. b, The residuals show that the contamination from the stellar 
photosphere is negligible in the observed spectrum. 
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for atypical noise level of ø = 0.15 mJy. The best-fitting emitting radius R for all emission at a temperature of about 600K. 
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Extended Data Fig. 5|Continuum-subtractedspectruminthe15ymregion used for illustrative purposes only. Rotational lines of H,O (J = 14515-1321, 
showing the detected Q-branch of CO, (orange). The shape of this feature is J=1469- 13310; Eu ~ 4300 K) are also detected and they are reasonably well 
sensitive to temperature and is well-fitted by an LTE slab model with 7 ~ 200K. reproduced by the best-fit model for the 7 um region (blue). This could indicate 
Thestrength of this feature can be reproduced with M(CO,) = 1.5 x 10” cm’ and that there is no additional reservoir of water at cooler temperature. 

R=0.1Au. However, the aforementioned parameters are degenerate and are 
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Extended Data Fig. 6| WISE Time-series photometry of PDS 70. Errorbars 
WISE 2 and WISE 4 indicating a’seesaw’like time variability (see Methods for 


represent 1s.d. a-b, WISE 3 (W3) is not anticorrelated with WISE land WISE 2 
due to the dominant 10 zm silicate emission which does not vary substantially further details). 


